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Abstract We present an algorithm of finding numerical 
solutions of pulsar equation. The problem of finding the 
solutions was reduced to finding expansion coefficients 
of the source term of the equation in a base of orthogo- 
nal functions defined on the unit interval by minimizing 
a multi-variable mismatch function defined on the light 
cylinder. We applied the algorithm to Scharlemann & 
Wagoner boundary conditions by which a smooth solu- 
tion is reconstructed that by construction passes success- 
fully the Gruzinov's test of the source function exponent. 

Keywords pulsars: general • stars: neutron • stars: 
rotation 
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1 Introduction 

Pulsar equation describes structure of electromagnetic 
fields and currents in magnetosphcrc of aligned rotator. 
The structure is uniquely determined by a scalar function 
W which is a solution of the equation 



(1 - p 2 )V 2 # - -d p W + Ftm = 0. 
p 



(1) 



The W and the unknown source function F(W) define 
structure of electromagnetic fields and currents in pulsar 
neighborhood. We use cylindrical coordinates p, z and 
cj) with the axis of rotation of the pulsar as the sym- 
metry axis. To solve the equation we assume Scharle- 
mann & Wagoner conditions which we specify later on. 
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This however means, that we do not consider the particle 
production issue nor we make any tests for the shape of 
emitting regions. For simplicity we do not take into ac- 
count return currents, however they can be incorporated 
by adding delta function in the source term. 



2 Assumptions of Scharlemann &: Wagoner 
model 

To describe structure of magnetosphere Scharlemann & 
Wagoner assumed in [6] the following conditions 

1. Axial symmetry and stationarity The assump- 
tions follow naturally from the simplifying require- 
ment that the axis of rotation of a neutron star over- 
laps with its dipole momentum (aligned rotator) . The 
axis is therefore the axis of cylindrical symmetry of 
the whole system and the system is independent of 
time. The model is a particular case of the oblique ro- 
tator ([5] and [7]). Consequently, the system does not 
generate the pulsar effect because there is no modula- 
tion. However, we consider the case of aligned rotator 
for simplicity. Such issue is interesting, too, and it was 
discussed in many articles (like [2], [6]). 

2. Force- free approximation In this approximation 
one assumes that inertial and gravitational forces are 
negligible, by which only electric and magnetic forces 
are taken into account. The assumption agrees with 
considerations of Goldreich & Julian [3] who showed 
that the electromagnetic interaction are stronger than 
gravitational forces by a factor of 10 8 for protons, 
and 10 11 for electrons. They also found that charge 
density just above the neutron star surface cannot be 
zero. One can imagine pulsar as extremely big atomic 
nucleus. 

3. Rigid co-rotation Although we are aware of the 
current works by Timokhin ([5] [8]) who considers 
models with critical point located inside the light 
cylinder, we present our algorithm assuming that the 
point is located just on the cylinder as in Scharle- 
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mann & Wagoner model. However, boundary condi- 
tions can always be adapted to account for Timokhin 
assumptions. Additionally, one may take into account 
return currents but we do not consider them in our 
short paper aimed at presenting our algorithm. 
Asymptotics It is assumed that asymptotically \P 
tends to the profile of the split magnetic monopole of 
unit charge, while in the vicinity of the center a mag- 
netic dipole should be located. The solution should be 
smooth everywhere apart from the equatorial plane 
outside the light cylinder. 



3 Numerics 

To find W with Scharlemann & Wagoner conditions in 
the whole physical space we adapted pulsar equation as 
follows 

1. Compactification is required to perform calcula- 
tions on a finite size lattice. For example, one can 
choose the mapping 

2p 2z 



y 



(2) 



+ 7 

that transforms the original infinite physical domain 
(p, z) £ [0, +oo] x [— oo, +00] onto the square (x, y) € 
[0,2] x [—2,2]. By symmetry it suffices to consider 
only the region y > 0. As an example in figure [1] 
there are shown profiles of the monopole and of the 
dipole in the compactified domain. 
Boundary conditions Scharlemann & Wagoner bound- 
ary conditions [6] in the compactified domain are 
given in table below 
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In fact, we could use any boundary conditions with 
our algorithm, in particular, the Y point could be 
moved to the interior of the light cylinder. 
Discretization Once we discretize our domain we 
may integrate the pulsar equation. For a given n the 
compactified lattice is defined by nodal points of the 
grid (j,k) such that x — 2j/{2n + 1), y = k/n, j 
0, 1, . . . , 2n + 1, k = 0, 1, . . .2n and U jtk = &(x, y), 
thus the singular line x = I is not used during inte- 
gration. We thereby avoided the cumbersome prob- 
lem of matching solutions along the light cylinder and 
carried on calculations on a single grid. The integra- 
tion grid together with boundary conditions is shown 
in figure O 



The source function The source function F(\P) is 
nonzero in the domain <P = (0, 1) (such normaliza- 
tion is possible) , then in regions where \P > 1 corota- 
tion takes place. Morover, from the analysis of pul- 
sar equations with Scharlemann & Wagoner condi- 
tions it follows that F'(0) — 4. Additionally, we take 
into account the result of paper by Gruzinov [4] that 
F(W) - (1 - <f) Q with a ps 7/12 as & -> 1". There- 
fore, one may expand F(<F) in the basis of the Jacobi 

polynomials fi(^) — P/^ 6 ' 4 (<f r ), i = Q, 1, 2, . . ., where 
we used the convention that F{fP) is the same as for 
the unit charge monopole at Cj = 0, i = 0, 1,2,.. .. 

F(W) = 2<?(1 -¥)(2-¥) + S l c l ^ 2 {\ - &) a fi{V){3) 

The key idea is to find the expansion coefficients such 
that be smooth on the light cylinder. This can be 
done by minimizing an error function £ (defined be- 
low) measuring departure from smoothness. By tak- 
ing only a few initial terms of the expansion the prob- 
lem of finding a solution is reduced to finding the 
minimum of £ which is standard in numerical anal- 
ysis. As an aside, we remark that by adding a con- 
tinuous representation of a delta function S n to the 
above definition of i ? (!f r ), that is, by replacing F(<F) 
with F(ip) + /3S n (l — <?"), on can similarly find a so- 
lution with return currents along separatrix, where /3 
is the additional parameter to be determined by the 
minimization. 

The error function To compute the error function 
we used the formula 



£(c-i,c , Ci, 



(4) 
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which utilizes the original smoothness condition 2d p l l / ■ 
F^) at p — 1, and Wk is arbitrary discrete weight 
function. The c_i is proportional to the dipole mo- 
mentum of \P and Co , ci , . . . are the expansion coeffi- 
cients. The key point of our algorithm is, for a given 
to, to find such a point {c_i, Co, c±, Cm,} at which 
£ (c-i, Co, ci, Cm) has the minimum. Finding min- 
ima of multi- variable functions is a standard problem 
in numerical analysis. It should be clear that the al- 
gorithm of finding solutions differs qualitatively from 
the method presented by CKF in [2]. 



= 4 The results 



In picture [2J and 0] are shown the whole profile of <F in 
the compactified domain and the source function, respec- 
tively, which were found with the help of our algorithm 
for n — 15 (we used small n, where n is the grid size, to 
show that quite good results can be obtained with the 
help of our algorithm even on very small size grids). In 



(5) 



An algorithm for solving the pulsar equation 



3 




Fig. 1 The monopole and the dipole in the compactified domain for y > 0. 




Fig. 2 The solution of pulsar equation with Scharlemann & Wagoner conditions shown in the compactified domain. 



the table below [4] there are shown also the correspond- 
ing parameters cq and c\, . . . , c@ obtained for n = 15. We 
remind we neglected in the presentation the singular re- 
turn currents to speed-up finding solutions, but one can 
easily modify the expansion of F^) to account for it as 
discussed earlier. 



improve the results qualitatively. We can repeat our cal- 
culations with Timokhin ([5]) boundary conditions. 

In our solution the surface of \P is smooth apart from 
the equatorial plane outside the light cylinder. 
The introduction to the problem can be found in pQ. 



5 Discussion and summary 

The new method of calculating source function seems to 
work well. For bigger grids we don't presume they would 
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Fig. 3 Compaetified integration grid which covers whole physical space. The light cylinder which is singular surface of the 
pulsar equation is located between adjacent grid columns. 
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0.0005 


4 
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Table 1 Location of the minimum of error function [4] we 
found for n = 15 and m — 6; d are coordinates of the mini- 
mum and o i are uncertainties of the coordinates. 
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Fig. 4 Source function F(SI/) found for n — 15 and m — 6. 



